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Ultra High Energy Cosmic Rays from Compact Sources 
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The clustering of ultra high energy (above 10 20 eV) cosmic rays (UHECR) suggests that they 
might be emitted by compact sources. Statistical analysis of Dubovsky et al. (Phys. Rev. Lett. 
85 (2000) 1154) estimated the source density. We extend their analysis to give also the confidence 
intervals for the source density using a.) no assumptions on the relationship between clustered and 
unclustered events; b.) nontrivial distributions for the source luminosities and energies; c.) the 
energy dependence of the propagation. We also determine the probability that a proton created at a 
distance r with energy E arrives at earth above a threshold E c . Using this function one can determine 
the observed spectrum just by one numerical integration for any injection spectrum. The observed 14 
UHECR events above 10 20 eV with one doublet gives for the source densities 180+ 2 gg° • 10" 3 Mpc" 3 
(on the 68% confidence level) . We present detailed results for future experiments with larger UHECR 
statistics. 



I. INTRODUCTION 

The interaction of protons with photons of the cosmic 
microwave background predicts a sharp drop in the cos- 
mic ray flux above the GZK cutoff around 5 ■ f0 19 eV 
0. The available data shows no such drop. About 20 
events above 1 20 eV were observed by a number of ex- 
periments such as AGASA Fly's Eye ||], Haverah 
Park f§, Yakutsk || and HiRes §]. Since above the 
GZK energy the attenuation length of particles is a few 
tens of megaparsecs if an ultra high energy cosmic 
ray (UHECR) is observed on earth it must be produced 
in our vicinity (except for UHECR scenarios based on 
weakly interacting particles, e.g. neutrinos). 

Usually it is assumed that at these high energies the 
galactic and extragalactic magnetic fields do not affect 
the orbit of the cosmic rays, thus they should point back 
to their origin within a few degrees, fn contrast to the low 
energy cosmic rays one can use UHECRs for point-source 
search astronomy. (For an extragalactic magnetic field of 
fiG rather than the usually assumed nG there is no di- 
rectional correlation with the source ||. Note, that due 
to the Local Supercluster the magnetic field is presum- 
ably not less than 10 nG which results in a Larmor radius 
of few tens of megaparsecs for protons above 10 20 eV.) 
Though there are some peculiar clustered events, which 
we discuss in detail, the overall distribution of UHECR 
on the sky is practically isotropic. This observation is 
rather surprising since in principle only a few astrophys- 
ical sites (e.g. active galactic nuclei || or the extended 
lobes of radio galaxies |Hj) are capable of accelerating 
such particles, nevertheless none H] of the UHECR evets 
came from these directions. Hence it is generally believed 
jl2| that there is no conventional astrophysical explana- 
tion for the observed UHECR spectrum. 



There are several ways to look for the source inhomo- 
geneity from the energy spectrum and spatial directions 
of UHECRs. One possibility is to assume that the source 
density of UHECRs is proportional to the galaxy densi- 
ties jL3j. Another approach is to analyze the clustering 
properties of the unknown sources by some correlation 
length @. 

Clearly, the arrival directions of the UHECRs mea- 
sured by experiments show some peculiar clustering: 
some events are grouped within ~ 3°, the typical an- 
gular resolution of an experiment. Above 4 • 10 19 eV 92 
cosmic ray events were detected, including 7 doublets and 
2 triplets. Above 10 20 eV one doublet out of 14 events 
were found p5| . The chance probability of such a clus- 
tering from uniform distribution is rather small fll5] , |l6| |. 
(Taking the average bin 3° the probability of generating 
one doublet out of 14 events is 11%.) 

The clustered features of the events initiated an in- 
teresting statistical analysis assuming compact UHECR 
sources |l7|]. The authors found a large number, ~ 400 for 
the number of sources^] inside a GZK sphere of 25 Mpc. 
They assumed that 

a.) the number of clustered events is much smaller than 
the total number of events (this is a reliable assump- 



* approximately 400 sources within the GZK sphere results in 
one doublet for 14 events. The order of magnitude of this re- 
sult is in some sense similar to that of a "high-school" exercise: 
what is the minimal size of a class for which the probability of 
having clustered birthdays -at least two pupils with the same 
birthdays- is larger than 50%. In this case the number of 
"sources" is the number of possible birthdays ~ 400. In order 
to get the answer one should solve 365!/ [365 fe (365 -fc)!] < 0.5, 
which gives as a minimal size k=23. 
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tion at present statistics; however, for any number of 
sources the increase of statistics, which will happen in 
the near future, results in more clustered events than 
unclustered) , 

b. ) all sources have the same luminosity which gives 
a delta function for their distribution (this unphysical 
choice represents an important limit, it gives the small- 
est source density for a given number of clustered and 
unclustered events) 

c. ) The GZK effect makes distant sources fainter; how- 
ever, this feature depends on the injected energy spec- 
trum and the attenuation lengths and elasticities of the 
propagating particles. In jL7 an exponential decay was 
used with an energy independent decay length of 25Mpc. 

In our approach none of these assumptions are used. 
In addition we include spherical astronomy corrections 
and in particular give the upper and lower bounds for the 
source density at a given confidence level. As we show the 
most probable value for the source density is really large; 
however, the statistical significance of this result is rather 
weak. At present the small number of UHECR events 
allows a 95% confidence interval for the source density 
which spreads over four orders of magnitude. Since future 
experiments, particularly Pierre Auger fll8| , |i9"| ] , will have 
a much higher statistical significance on clustering (the 
expected number of events of 10 20 eV and above is 60 
per year 20 1), we present our results on the density of 
sources also for larger number of UHECRs above 10 20 
eV. 

In order to avoid the assumptions of Jl7| a combined 
analytical and Monte-Carlo technique will be presented 
adopting the conventional picture of protons as the ul- 
tra high energy cosmic rays. Our analytical approach of 
Section gives the event clustering probabilities for any 
space, luminosity and energy distribution of the sources 
by using a single additional function P(r, E; E c ), the 
probability that a proton created at a distance r with 
energy E arrives at earth above the threshold energy E c 

we 



With our Monte-Carlo technique of Section [II 



determine the probability function P(r, E; E c ) for a wide 
range of parameters. Our results for the present and fu- 
ture UHECR statistics are presented in Section IV We 
summarize in Section [v|. 



II. ANALYTICAL APPROACH 

The key quantity for finding the distribution functions 
for the source density is the probability of detecting k 
events from one randomly placed source. The number of 
UHECRs emitted by a source of A luminosity during a 
period T follows the Poisson distribution. However, not 
all emitted UHECRs will be detected. They might loose 
their energy during propagation or can simply go to the 
wrong direction. 

For UHECRs the energy loss is dominated by the 
pion production in interaction with the cosmic microwave 



background radiation. In rcf. the probability func- 
tion P(r, E, E c ) was presented for three specific threshold 
energies. This function gives the probability that a pro- 
ton created at a given distance from earth (r) with some 
energy (E) is detected at earth above some energy thresh- 
old (E c ). The resulting probability distribution can be 
approximated over the energy range of interest by a func- 
tion of the form 

P(r, E, E c ) « exp[-a(£ c )r 2 ex V (b(E c ) / E)} (1) 

The appropriate values of a and b for i? c /(10 20 eV) =1,3, 
and 6 are, respectively a/(l(U 4 Mpc~ 2 ) =1.4, 9.2 and 11, 
6/(10 20 eU) =2.4, 12 and 28. 

For the sources we use the second equatorial coordinate 
system: x is the position vector of the source character- 
ized by (r, S, a) with 5 and a beeing the declination and 
right ascension, respectively. The features of the Poisson 
distribution enforce us to take into account the fact that 
the sky is not isotropically observed. There is a circum- 
polar cone, in which the sources can always be seen, with 
half opening angle 6' (5' is the declination of the detec- 
tor, for the experiments we study 6' ps 40° — 50°). There 
is also an invisible region with the same opening angle. 
Between them there is a region for which the time frac- 
tion of visibility, j(S,8') is a function of the declination 
of the source. It is straightforward to determine j(5, 5') 
for any 5 and S': 



j(5,5') = 




-tt/2 < 5 < 5' -tt/2 
(tan 8' tan 8) /it 
8' - tt/2 < 8 < tt/2 - 5' 
tt/2 -8' < 8 < tt/2 



(2) 



To determine the probability that a particle arriving from 
random direction at a random time is detected we have 
to multiply j(8, 8') by the cosine of the zenith angle 9. In 
the following we will use the time average of this function: 



V(S, 8') 



1 

T 



j(S, 8') -cos6>(<5, S',t)dt 



(3) 



Since 8' is constant, in the rest of the paper we do not 
indicate the dependence on it. Neglecting these spherical 
astronomy effects means more than a factor of two for the 
prediction of the source density. 

The probability of detecting k events from a source at 
distance r with energy E can be obtained by including 
P(r, E , E c )Arj(5) / (4-iTr 2 ) in the Poisson distribution: 



, „ .v exp[-P(r,E,E c )r)(8)j/r 2 ] 

[P(r,E,E c ) V (S)j/r 2 ]\ 



(4) 



where we introduced j — XT A/ (Air) and Ar](8) / {Airr 2 ) 
is the probability that an emitted UHECR points to a 
detector of area A. We denote the space, energy and 
luminosity distributions of the sources by p(x) , c{E) and 
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h(j), respectively. The probability of detecting k events 
above the threshold E c from a single source randomly 
positioned within a sphere of radius R is 

/> />oo />oo 

P k = / dV p(x) / dE c{E) / dj h(j) x 
Js R Je c Jo 

exp[-P(r,E,E c )r)(5)j/ 



k\ 



'"' [P(r,E,E c ) v (S)j/r 2 ] k . (5) 



Denote the total number of sources within the sphere 
of sufficiently large radius (e.g. several times the GZK 
radius) by N and the number of sources that gave k de- 
tected events by Nk- Clearly, N = J2^N and the total 
number of detected events is N e = J^o The proba- 

bility that for N sources the number of different detected 
multiplets are Nk is: 



P(N,{N k })=N!fl 1 Lp»>«. 



k=0 



(6) 



The value of P(N, {Nk}) is the most important quantity 
in our analysis of UHECR clustering. For a given set 
of unclustered and clustered events (iVi and N 2 ,N 3 ,...) 
inverting the P(N, {Nk}) distribution function gives the 
most probable value for the number of sources and also 
the confidence interval for it. If we want to determine the 
density of sources we can take the limit R — > oo, N — > oo, 
while the density of sources S = N/(^R 3 n) is constant. 

In order to illustrate the dominant length scale it is 
instructive to study the integrand /fc(r) of the distance 
integration in eqn. (||) 



Pi- 



„ J'*" - 



/poo poo 
dfip(x) / dE c(E) / dj h(j) 
JE a Jo 

CX p[-P(r,E,E c ) V (6)j/ 



k\ 



'"' [P(r,E,E c ) V (S)j/r 2 ] k . (7) 



Fig. [T] shows that /i(r), which leads to singlet events, is 
dominated by the distance scale of 10-15 Mpc, whereas 
/2(r), which gives doublet events, is dominated by the 
distance scale of 4-6 Mpc|]. These typical distances partly 



t It is interesting that the dominant distance scale for singlet 
events is by an order of magnitude smaller than the attenua- 
tion length of the protons at these energies (l a ~ 110 Mpc). 
This surprising result can be illustrated using a simple approx- 
imation. Assuming that the probability of detecting a parti- 
cle coming from distance r is proportional to exp(—r/l a )/r 2 , 
Pi will be proportional to J dfldrr 2 ■ exp[— j exp(— r/l a )/r 2 ] ■ 
exp(—r/l a )/r 2 . For the typical j values the r integrand has a 
maximum around 15 Mpc and not at l a - 



justify our assumption of neglecting magnetic fields. The 
deflection of singlet events due to magnetic field does not 
change the number of multiplets, thus our conclusions. 
The typical distance for higher multiplets is quite small, 
therefore deflection can be practically neglected. Clearly, 
the fact that multiplets are coming from our "close" 
neighbourhood does not mean that the experiments re- 
flect just the densities of these distances. The overwhelm- 
ing number of events are singlets and they come from 
much larger distances. Note, that these /i(r) and /bO") 
functions are obtained with our optimal j* value (cf. Fig. 
^ and explanation there and in the corresponding text). 
Using the largest possible j* value allowed by the 95% 
confidence region the dominant distance scales for fi(r) 
and /2(f) functions turn out to be 30 Mpc and 20 Mpc, 
respectively. 

Note, that Pk and then P(N, {Nk}) are easily deter- 
mined by a well behaving four-dimensional numerical in- 
tegration (the a integral can be factorized) for any c(E), 
h(j) and p(r) distribution functions . In order to illus- 
trate the uncertainties and sensitivities of the results we 
used a few different choices for these distribution func- 
tions. 

For c(E) wc studied three possibilities. The most 
straightforward choice is the extrapolation of the 'con- 
ventional high energy component' oc E~ 2 . Another pos- 
sibility is to use a stronger fall-off of the spectrum at 
energies just below the GZK cutoff, e.g. oc E~ 3 . These 
choices span the range usually considered in the litera- 
ture and we will study both of them. The third possibility 
is to assume that topological defects generate UHECRs 
through production of superheavy particles q. Accord- 
ing to (2j| these superheavy particles decay into quarks 
and gluons which initiate multi-hadron cascades through 
gluon bremstrahlung. These finally hadronize to yield 
jets. The energy spectrum, first calculated in for the 
Standard Model and in |2j| for the Minimal Supersym- 
metric Standard Model, in this case can be estimated by 
the function obtained from the HERWIG QCD genera- 
tor: 



c{x) = Ci 



exp[c 2v /ln(l/ a: )](l-a;) 2 
x^J\n(l/x) 



(8) 



where x = E/mx the ratio of the e nerg y and the mass 
of the decaying particle. The best fit |25| to the observed 
UHECR spectrum gives mx ~ 10 12 GeV for the mass of 
the decaying particle. This corresponds to ci w 0.0086 
and C2 ~ 2.77. We will use these ci and C2 values for our 
third choice of energy distribution, c(E). 



*Note, that these particles are not superheavy DM particles 
[^3), which are located most likely in the halo of our galaxy. 
These superheavy DM particles can also be considered as pos- 
sible sources of UHECR p^-pr| with anisotrpies in the arrival 
direction l|27f. 
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10 15 
r [Mpc] 

FIG. 1. The distributions fi(r) -solid line- and /2(f) 
-dashed line- of eqn. (j^). The singlet and doublet events 
are dominated by distance scale of 10-15 Mpc and 3-5 Mpc, 
respectively. 



In ref. [[l7| the authors have shown that for a fixed set of 
multiplets the minimal density of sources can be obtained 
by assuming a delta- function distribution for h(j). We 
studied both this limiting case (h(j) = S(j — j'*)) and a 
more realistic one with Schechter's luminosity function 



h(j)dj = h ■ (j/j* 



-1.25 



exp(-j/j*)<2(j/j*)- 



(9) 



The space distribution of sources can be given based 
on some particular survey of the distribution of nearby 
galaxies |l3| ] or on a correlation length ro characterizing 
the clustering features of sources . For simplicity the 
present analysis deals with a homogeneous distribution 
of sources randomly scattered in the universe (Note, that 
due to the Local Supercluster the isotropic distribution 
is just an approximation.). 

Fig. H shows the resulting Pk(j*) probability functions 
for the different choices of c(E) and h(j). The overall 
shapes of them are rather similar; nevertheless, relatively 
small differences lead to quite different predictions for the 
UHECR source density. The "shoulders" of the curves 
with Dirac-delta luminosity distributions got smoother 
for the Schechter's distribution. The scales on the fig- 
ures are chosen to cover the 98% confidence regions (see 
section IV for details). 

Note, that - assuming that UHECRs point back to 
their sources - our clustering technique discussed above 
applies to practically any models of UHECR (e.g. neutri- 
nos). One only needs a change in the P(r, E, E c ) proba- 
bility distribution function (e.g. neutrinos penetrate the 
microwave background uninhibited) and use the h(j) and 
c(E) distribution function of the specific model. 



10 1 10 
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1 10 1 io : 



J. [(Mpc) 2 



FIG. 2. The individual Pk(j*) functions for the different 
c(E) and h(j) choices. The column on the left corresponds to 
the Dirac-delta distribution h(j) = 5(j — j*), whereas the col- 
umn on the right shows the results for Schechter's luminosity 
distribution. The Erst, second and third rows correspond to 
the c(E) functions proportional to E~ 2 , E~ s and the super- 
heavy decay mode, respectively (see text). On each panel the 
individual lines from top to bottom are: 1 — Po, Pi, P2, Pz, 
P4, and P 5 . 
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III. MONTE-CARLO STUDY OF THE 
PROPAGATION 

Our Monte-Carlo model of UHECR studies the propa- 
gation of UHECR. The analysis of pll showed that both 
AGASA and Fly's Eye data demonstrated a change of 
composition, a shift from heavy -iron- at 10 17 eV to light 
-proton- at 10 19 eV Thus, the chemical composition of 
UHECRs is most likely to be dominated by protons. In 
our analysis we use exclusively protons as UHECR par- 
ticles, (for suggestions that air showers above the GZK 
cutoff are induced by neutrinos see p3| -) 

Using the pion production as the dominant effect of 
energy loss for protons at energies > 10 19 eV ref. Q] cal- 
culated P(r, E, E c ), the probability that a proton created 
at a given distance (r) with some energy (E) is detected 
at earth above some energy threshold (E c ). For three 
threshold energies the authors of [Q gave an approxi- 
mate formula, which we used in the previous section. 

In our Monte-Carlo approach we determined the prop- 
agation of UHECR on an event by event basis. Since 
the inelasticity of Bethe-Heitler pair production is rather 
small (« 10 -3 ) we used a continuous energy loss ap- 
proximation for this process. The inelasticity of pion- 
photoproduction is much higher (« 0.2 — 0.5) in the en- 
ergy range of interest, thus there are only a few tens of 
such interactions during the propagation. Due to the 
Poisson statistics of the number of interactions and the 
spread of the inelasticity, we will see a spread in the en- 
ergy spectrum even if the injected spectrum is mono- 
energetic. 

In our simulation protons are propagated in small steps 
(10 kpc), and after each step the energy losses due to pair 
production, pion production and the adiabatic expansion 
are calculated. During the simulation we keep track of 
the current energy of the proton and its total displace- 
ment. This one avoids performing new simulations for 
different initial energies and distances. The propagation 
is completed when the energy of the proton goes below a 
given cutoff. For the proton interaction lengths and in- 
elasticities we used the values of pl| , ^2[ . The deflection 
due to magnetic field is not taken into account, because 
it is small for our typical distances illustrated in Fig. [I]. 
This fact justifies our assumption that UHECRs point 
back to their sources (fo r a recent Monte-Carlo analysis 
on deflection see e.g. |p3[ ). 

Since it is rather practical to use the P(r, E, E c ) proba- 
bility distribution function we extended the results of |TJ] 
by using our Monte-Carlo technique for UHECR propa- 
gation. In order to cover a much broader energy range 
than the parametrization of ([!]) we used the following 
type of function 



P(r, E, E c ) = exp [-a ■ (r/1 Mpc) fc 



(10) 



Fig. U demonstrates the reliability of this parametriza- 
tion. The direct Monte-Carlo points and the fitted func- 
tion (eqn. @ with a = 0.0019 and b = 1.695) are 




40 60 80 
r [Mpc] 

FIG. 3. The direct Monte-Carlo points and the fitted func- 
tion P{r,E,E c ) = exp [-a-(r/ IMpc) 6 ] for E c = 10 20 eV 
and E = 2 ■ 10 20 eV. The fitted curve corresponds to 
a = 0.0019 and b = 1.695. 



plotted for E c = 10 20 cV and E = 2 • 10 20 cV. Fig. | 
shows the functions a(E/E c ) and b(E/E c ) for a range of 
three orders of magnitude and for five different thresh- 
old energies. Just using the functions of a(E/E c ) and 
b(E/E c ), thus a parametrization of P(r, E, E c ) one can 
obtain the observed energy spectrum for any injection 
spectrum without additional Monte-Carlo simulation. 



IV. RESULTS 



In order to determine the confidence intervals for the 



source densities we used the frequentist method |34| . We 
wish to set limits on S, the source density. Using our 
Monte-Carlo based P(r, E, E c ) functions and our ana- 
lytical technique we determined p(Nx, AT 2 , Ns, S; j*), 
which gives the probability of observing iVi singlet, N2 
doublet, A 3 triplet etc. events if the true value of the den- 
sity is S and the central value of luminosity is j* . The 
probability distribution is not symmetric and far from 
being Gaussian. For a given set of {Ni,i = 1,2, ...} the 
above probability distribution as a function of S and j* 
determines the 68% and 95% confidence level regions in 
the S — j* plane. Fig. || shows these regions for our "fa- 
vorite" choice of model (c(E) cx E~ 3 and Schechter's lu- 
minosity distribution) and for the present statistics (one 
doublet out of 14 UHECR events). The regions are 
deformed, thin ellipse-like objects in the log(j*) versus 
log(S') plane. Since is a completely unknown and in- 
dependent physical quantity the source density can be 
anything between the upper and lower parts of the con- 
fidence level regions. For this model our final answer 



for the density is 180 



+2730(8817) 
-165(174) 



10 3 Mpc 3 , where 
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FIG. 4. The functions a(E/E c ) -left panel- and b(E/E c ) 
-right panel- for the probability distribution function 
P(r, E, E c ) using the parametrization exp[— a-(r/l Mpc) b ] for 
Eve different threshold energies (5-10 19 eV, 10 20 eV, 2-10 20 eV, 
5 • 10 20 eV and 10 21 eV). 




10 1 10* 
J. [(Mpc) 21 

FIG. 5. The la (68%) and 2a (95%) confidence level 
regions for j* and the source density (14 UHECR with one 
doublet). The most probable value is represented by the tri- 
angle. The upper and lower boundaries of these regions give 
for the source density 180^ 2 gg^^ 7 ' 
(95% )confidence level. 



10" 3 Mpc~ 3 on the 68% 



the first errors indicate the 68%, the second ones in the 
parenthesis the 95% confidence levels, respectively. The 
choice of [|l7) -Dirac-delta like luminosity distribution- 
and, for instance, conventional E~ 2 energy distribution 
gives much smaller value: 2.77+f 5 3[ 2 ' 1 7 6 ( ] ) lCr 3 Mpc~ 3 . For 
other choices of c(E) and h(j) see Table |[ Our results for 
the Dirac-delta luminosity distribution are in agreement 
with the result of JlTj within the error bars. Neverthless, 
there is a very important message. The confidence level 
intervals are so large, that on the 95% confidence level 
two orders of magnitude smaller densities than suggested 
as a lower bound by |L7| are also possible. 

As it can be seen there is a strong correlation between 
the luminosity and the source density. Physically it is 
easy to understand the picture. For a smaller source 
density the luminosities should be larger to give the same 
number of events. However it is not possible to produce 
the same multiplicity structure with arbitrary luminosi- 
ties. Very small luminosities can not give multiplets at 
all, very large luminosities tend to give more than one 
doublet. 

The same technique can be applied for any hypotheti- 
cal experimental result. For fixed {Nk} the above proba- 
bility function determines the 68% confidence regions in 
S and j*. Using these regions one can tell the 68% con- 
fidence interval for S. The most probable values of the 
source densities for fixed number of multiplets are plot- 
ted on Fig. ^| with the lower and upper bounds. The 
total number of events is shown on the horizontal axis, 
whereas the number of multiplets label the lines. Here 
again, our " favorite" choice of distribution functions were 
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used: c(E) oc E~ 3 and h(j) of eqn. @. 

It is of particular interest to analyze in detail the 
present experimental situation having one doublet out of 
14 events. Since there are some new unpublished events, 
too, we study the a hypothetical case of one or two dou- 
blets out of 24 events. The 68% and 95% confidence 
level results are summarized in Table | for our three en- 
ergy and two luminosity distributions. It can be seen 
that Dirac-delta type luminosity distribution really gives 
smaller source densities than broad luminosity distribu- 
tion, as it was proven by Jl7| . Less pronounced is the 
effect on the energy distribution of the emitted UHE- 
CRs. The c(E) oc E~ 3 case gives somewhat larger values 
than the other two choices (c(E) oc E~ 2 or given by the 
decay of a superheavy particle) . The confidence intervals 
are typically very large, on the 95% level they span 4 or- 
ders of magnitude. An interesting feature of the results 
is that "doubling" the present statistics with the same 
clustering features (in the case studied by the table this 
means one new doublet out of 10 new events) reduces the 
confidence level intervals by an order of magnitude. The 
reduction is far less significant if we add singlet events 
only. Inspection of Fig. |] leads to the coclusion that 
experiments in the near future with approximately 200 
UHECR events can tell at least the order of magnitude 
of the source density. 

V. SUMMARY 

We presented a technique in order to statistically an- 
alyze the clustering features of UHECR. The technique 
can be applied for any model of UHECR assuming small 
deflection. The key role of the analysis is played by the 
Pk functions defined by eqn. (|^), which is the probability 
of detecting k events above the threshold from a single 
source. Using a combinatorial expression of eqn. (g) the 
probability distribution for any set of multiplets can be 
given as a function of the source density. 

We discussed several types of energy and luminosity 
distributions for the sources and gave the most probable 
source densities with their confidence intervals for present 
and future experiments. 

The probability P(r, E, E c ) that a proton created at a 
distance r with energy E arrives above the threshold E c 
is determined and parametrized for a wide range of 
threshold energies. This result can be used to obtain the 
observed energy spectrum of the UHECR for arbitrary 
injection spectrum. 

In ref. jlTj the authors analyzed the statistical fea- 
tures of clustering of UHECR, which provided constraints 
on astrophysical models of UHECR when the number of 
clusters is small, by giving a bound from below. In our 
paper we have shown that there is some constraint, but it 
is far from being tight. At present statistics the 95% con- 
fidence level regions usually span 4 orders of magnitude. 
Two orders of magnitude smaller numbers than the pre- 
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FIG. 6. The most probable values for the density of 
sources as a function of the total number of events (mid- 
dle panel). The number of multiplets are indicated on the 
individual lines in the form: N%, N3, N4, where A^A^ and 
N4 represent the appropriate values for doublets, triplets and 
quartets. The upper and lower panels correspond to the 84 
percentile and 16 percentile lines (upper and lower bounds of 
the 68% confidence intervals), respectively. 
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TABLE I. The most probable values for the source densi- 
ties and their error bars given by the 68% and 95% conBdence 
level regions (the latter in parenthesis). The numbers are in 
units of 10 -3 Mpc -3 The three possible energy spectrums are 
given by a distribution proportional to E~ 2 , E~ 3 , or by the 
decay of a 10 12 GeV particle (denoted by "decay"). The lu- 
minosity distribution can be proportional to a Dirac-delta or 
to Schechter's luminosity function (denoted by "SLF"). Re- 
sults are listed for the observed 1 doublet out of 14 events and 
for two hypothetical cases (1 doublet out of 24 events and 2 
doublets out of 24 events). 



diction of jlTj (their eqn. (13) suggests for the density of 
sources ~ 6 ■ 10 -3 Mpc -3 ) can also be obtained. Adding 
10 new events with an additional doublet the confidence 
interval can be reduced to 3 orders of magnitude and the 
increase of the UHECR events to 200 can tell at least the 
order of magnitude of the source density. 
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